home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
C/C++ Users Group Library 1996 July
/
C-C++ Users Group Library July 1996.iso
/
listings
/
v_09_05
/
9n05050a
< prev
next >
Wrap
Text File
|
1991-02-05
|
3KB
|
157 lines
/* fast Fourier transform (FFT)
from Handbook of C tools for Scientists and Engineers by L. Baker
DEPENDENCIES:
header file complex.h required
*/
#include <stdio.h>
#include "complex.h"
#define twopi 6.283185307
int iterp;/* global to return count used*/
#define max(a,b) (((a)>(b))? (a): (b))
#define abs(x) ((x)? (x):-(x))
#define DOFOR(i,to) for(i=0;i<to;i++)
main(argc,argv) int argc;char **argv;
{int i,ii,nh,n=16;
double invn,exp(),dt=.25,omega,realpt,imagpt;
struct complex w[32],wi[32], data[32];
nh=n>>1;
DOFOR(i,n)
{
CMPLX(data[i], exp(-dt*(i)),0.);
};
/* caveat*/ data[0].x=.5;/*not 1. see text*/
printf(" before transform:\n");
DOFOR(i,n){printc(&(data[i])) ;printf("\n");}
invn=1./n;
fftinit(w,wi,n);
printf(" w factors:\n");
DOFOR(i,n){printc(&(w[i])) ;printc(&(wi[i])) ;printf("\n");}
printf(" transformed:\n");
fft(data,w,n);
DOFOR(i,n){printc(&(data[i])) ;printf("\n");}
printf(" scaled by dt and compared to analytic answer:\n");
DOFOR(i,n)
{
ii= (i-nh) ? i : n-i ;
omega= twopi*ii/(n*dt);
realpt= 1./(1.+omega*omega);
imagpt= -realpt*omega;
printf(" %f %f %f %f\n", dt*data[i].x,dt*data[i].y,realpt,imagpt);
}
/* inverse fft*/
fft(data,wi,n);
printf(" transformed back and scaled by 1/N:\n");
DOFOR(i,n)
{
CTREAL(data[i],data[i],(invn));
printc(&(data[i]));printf("\n");
}
exit(0);
}
int bitr(k,logn) int k,logn;
{
int ans,j,i;
ans=0;
j=k;
DOFOR(i,logn)
{
ans=(ans<<1)+(j&1);
j=j>>1;
}
return(ans);
}
int log2(n) int n;
{
int i;
i=-1;/* will return -1 if n<=0 */
while(1)
{
if(n==0)break;
n=n>>1;
i++;
}
return(i);
}
printc(x) struct complex *x;
{
printf("%f %f",x->x,x->y);
return;
}
fftinit(w,wi,n) int n; struct complex w[],wi[];
{
int i;
double realpt,imagpt,cos(),sin();
double factr,angle;
factr=twopi/n;
DOFOR(i,n)
{angle=i*factr;
realpt=cos(angle);imagpt=sin(angle);
CMPLX(w[i],realpt,imagpt);
CMPLX(wi[i],realpt,(-imagpt));
}
return;
}
fft(x,w,n) int n; struct complex w[],x[];
{
int n1,logn,i,j,k,l,logl,p;
struct complex s,t;
logn=log2(n);
n1=n>>1;
j=logn-1;
/* transform*/
k=0;
DOFOR(logl,logn)
{
do{
DOFOR(i,n1)
{
p=bitr((k>>j),logn);
l=k+n1;
CONJG(s,w[p]);
CMULT(t,s,x[l]);
CSUB(x[l],x[k],t);
CADD(x[k],t,x[k]);
k++;
};/* dofor i*/
k+=n1;
}while (k<n);
k=0;
j--;
n1=n1>>1;
}
/*reorder*/
for(i=1;i<n;i++)
{
k=bitr(i,logn);
if (i>k)
{
/*exchange i,k elements*/
CLET(s,x[i]);
CLET(x[i],x[k]);
CLET(x[k],s);
}
};
return;
}